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ABSTRACT. The q-model of Coppersmith et al. has renewed interest in understanding the 
forces generated along the walls and at the bottom of a silo filled with a granular material. Fluc- 
tuations in the mean stress have been characterized for the q-model, and related to experimental 
work on stress chains. The classical engineering approach to bin loads follows from Janssen's 
analysis, predicting a saturation of stress, as a function of depth, in a tall silo. In this note we 
re-examine the Janssen theory, introducing randomness into the important parameters in the 
theory. The Janssen analysis relies on assumptions not met in practice. For this reason, we 
numerically solve the PDEs governing the equilibrium of forces in a bin, again including ran- 
domness in parameters. We show that the most important of these parameters is a coefficient of 
friction at the wall of the bin. This random friction model combines some features of fiuctuations 
as seen in experiments, with a classical continuum mechanics approach to describing granular 
materials. 
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1. INTRODUCTION The classical engineering theory of Janssen ^ provides an estimate 
for the mean vertical stress in a silo filled with a granular material. The principal feature of the 
Janssen analysis is that, under passive stress conditions, the mean stress saturates, asymptoting 
to a value depending on bin radius and wall and internal friction coefficients, but independent 
of the height. The Janssen theory relies on two assumptions, assumptions which do not hold 
in practice. Nevertheless, the analysis gives a reasonable estimate of the bin loads, and its 
simplicity is its virtue. Several analyses have attempted to remove some of the assumptions of 
the Janssen theory; the interested reader should consult [Q. 

Recently Coppersmith et.al. @, |] have developed a model for the force distribution in a bin. 
In this "q-model" , any one particle within the sample transmits its weight to neighbors that are 
below, in a random manner. The authors derive a mean field theory based on this model, and 
find fiuctuations in the forces felt by the lowest row of particles. Under most assumptions on 
the choice of random number distribution, the number of occurrences of a fiuctuation of a given 
size decays exponentially with size. 

The current incarnation of the q-model is scalar: only the vertical force is balanced. Recently, 
Socolar has introduced a generalization, the so-called a-model, which balances vertical and 
horizontal forces and angular moment. The a- model contains three random variables, and 
analysis appears difficult. However numerical simulations modeling rough walled bins appear 
consistent with the classical continuum theory; numerical simulations modeling infinitely wide 
bins appear consistent with the q-model for a special distribution of the q's. 

To provide a framework for introducing fiuctuations into a continuum setting, we incorporate 
some of the randomness of the q-model into the Janssen analysis. In Section 2 we reconsider the 
Janssen derivation and include a random component into the grain friction, and reformulate the 
balance law as a stochastic differential equation. Standard results of stochastic calculus provide 
an estimate of the mean stress and its variance, at any height. In Section 3, we numerically solve 
the complete stress equilibrium equations, assuming a Mohr-Coulomb constitutive relation, and 
again including a random component in the friction. Under passive loading, the stress saturates; 
stress fluctuations are not signiflcant until near saturation. 

An experimental flnding closely related to the current note is [|] . That paper reports careful 
measurements of force fluctuations in tall narrow bins, bins whose widths range from three 
to eight grain diameters and whose depth ranges up to about 100 grains diameters. Measured 
average vertical stress at any depth is systematically higher than predicted by the Jansen theory, 
and fluctuations in this stress range up to about 20%. These fluctuations are apparent only after 
the stress starts to saturate. (N.B. Socolar |jlO| also flnds the Janssen stress smaller than his 
calculated average stress, at any depth.) These experiments also demonstrated a dependence of 
stress on ambient temperature, an effect we do not consider here. 



2. GENERALIZED JANSEN ANALYSIS We briefly review Janssen's theory, and provide 
a stochastic generalization of that analysis. See |^ for the fundamental mechanics of granular 
media. All of this study is restricted to two space dimensions. 

Let the average vertical stress be denoted a = /f^//2 o''^^{x,y) dx, where a^^^a^^^a^^ are 
the XX, xy, yy-components, respectively, of the (symmetric) stress tensor T. 

Consider the force diagram in Figure 1; at equilibrium, the average stress at y and y + Ay, 
gravity, and wall friction f are balanced: 

dy<y + = pg ■ (1) 
Now we make two assumptions, critical to the Janssen theory, but which do not hold in practice: 
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Figure 1: Force balance for Janssen's analysis. On the slice, stresses and gravity are balanced 
by wall friction. 



1. At every point cr^'^ and a^^ are the principal stresses (i.e., the eigenvalues of the stress 
tensor) and the Coulomb frictional condition implies a^^{x,y) = K(7^y{x,y), K = j^, 
s = sin(0) and (p is the internal friction angle; 

2. Along the wall, f = a^^(±L)/2, y) = 5(t^^(±D/2, y) where 5 = tan(0^), 0^ is the wall- 
material friction angle 

Combining these assumptions, we arrive at the equation 

dya + aa = pg, a = . (2) 

Solving subject to a ^ as ?/ — > gives 

^{y) = ~ ~ exp{-ay)) . (3) 

It is apparent that the average stress saturates, the asymptotic value — depending on the 
material and wall parameters and the bin diameter. 

The formula for K is based on the assumption that the stress field is in the passive state, 
with the a;a;-stress the major principal stress (the larger of the eigenvalues) and the yy-stiess 
the minor (the smaller eigenvalue). If the material is in the active state, the yy-stress is major, 
the xx-stress minor, and K is replaced by K~^. For a typical material, may be 30°, so = 3 
in the passive state. In the active state this parameter is |, and saturation of the stress requires 
a bin that is an order of magnitude taller. 

Now assume the coefficient of the stress in (^ has a mean and a fluctuating component. 
This fluctuating component might arise from randomness in the friction angle, for example. 
Assuming an Ito formulation for the resulting stochastic differential equation, write 

da = — aady — eadW + pgdy . (4) 

Here dW{y) is a Wiener measure associated with the random fluctuations, and e is a measure 
of the size of the fluctuations. Standard arguments give the following results (see, e.g., chapt. 
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8 in A formal solution may be obtained by a variation of parameters argument, but more 
insightful are formulae for the first and second moments. The mean of the solution, m = S{a), 
is, not surprisingly, the Janssen solution (^. The second moment P = S{d''^), satisfies 

P = {-2a + e^)P + 2mpg . 

Thus 



i^ = ^^[l-exp(-(2a-e»] (5) 



The standard deviation is \/P — m^, and an order of magnitude estimate gives the deviation 
~ after the stress has saturated. 

An alternative hypothesis is that randomness in packing leads to fluctuations in the density, 
and thus to fluctuations in the stress. That is, the weight pg must include a random component 
due to voids. This assumption leads to the equation 

da = — aady + pgdy + epgdW . (6) 

The mean of the solution is, again, given by (H). The standard deviation is ^^^[1— exp(— 2q;?/)]^. 



3. EQUILIBRIUM ANALYSIS The Janssen analysis relies on assumptions not met in 
practice. In this section we solve the full stress equilibrium equations for a Coulomb material 
in a bin. Although analysis is possible in the limiting case of smooth walls (see [§), this section 
determines solutions numerically. 
Stress equilibrium is written 

dxf^xx + dya^y = (7) 
dxCTyx + dyayy = pg . (8) 

A common constitutive assumption is that the material is Mohr-Coulomb, at incipient yield. 
That is, one assumes the ratio of the shear stress, r, to the mean stress, a, is a constant, where 

a = r = (9) 

and ai, a2 are the eigenvalues of the stress tensor T. The Mohr-Coulomb condition reads 

l^s. (10) 

The Mohr-Coulomb condition can be viewed as a nonlinear relation for, say, ayy in terms of 
a^x and a^y It is often convenient to make a change of variables that incorporates this relation. 
With the mean stress a defined above, introduce the angle ip, measured from the horizontal, 
such that {cos{ip), sin('?/')) is an eigenvector of T associated with ai. Then write 



T 




cos(2'?/') sin(2-?/') 
sm{2ip) — cos(2?/') 



This equation specifies the stresses in terms of two dependent variables, a, tp, whose evolution 
is determined by the equihbrium equations. 
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This change of variables may be used to rewrite the momentum equations as 



1 + s cos(2-?/') —2as sm{2ip) 

ssm{2ip) 2as cos{2ip) 

f ssin(2-?/') 2oscos{2il)) 

\ 1 — scos(2^/') 2crs sin(2?/') 



y 





We nondimensionahze by scahng length by the bin diameter Z), and stress by pgD. All calcu- 
lations are reported in non-dimensional units. The independent variable are —1/2 < x < 1/2 
and < y < H. This system of PDEs is strictly hyperbolic, with characteristics inclined at 
an angle ±(| — |) from the direction of major principle stress. The y-direction may be taken 
as the time-like direction. "Initial" conditions for a and ip are imposed at the top of the fill, 
y = 0, and the equations solved downward. At the boundaries, the bin walls x = ±1/2, the 
wall friction angle is imposed: ip = 6. The system of equations is solved by a modification of 
the TVD/Central Difference scheme of Nessyahu and Tadmor [|lT|. The method is second-order 



accurate, and designed to avoid spurious oscillations common to many higher-order schemes for 
hyperbolic systems. For the computations reported here, a gridsize of Ax = 0.02 was used. 
On very coarse grids, fluctuations are larger than shown; after sufficient refinement, the size of 
fiuctuations appears to stabilize. 

To introduce fiuctuations, at each gridpoint at each level, the friction angle is chosen with 
a random component. Specifically, if is the nominal friction angle, the angle used is (pjiuct = 
0(1.0 + CO "where ^ is chosen randomly from a uniform distribution between [-0.5,0.5] and ( is 
an adjustable parameter measuring the extent of variation in the friction angle. Depending on 
testing apparatus, variations in measurements of the internal friction angle are as large as ±5°, 
more than 10% of typical values 0; without good measurements of the span of friction angles, 
we conservatively set ( = 0.1. So, for a nominal friction angle of 30°, (pjiuct G [28.5, 31.5]; the 
sine of this angle is used in the constitutive relation, and this sine G [.477, .522], about a ±5% 
swing. Of course variations of friction in a real sample may have spatial correlations; absent 
good modeling justification for a particular choice of correlation, none is used here. A random 
component of the wall friction angle (the boundary condition 6) is added in a manner similar 
to 0. We emphasize that our choice of ( sets the imposed variation in friction angle, and thus 
of the stress, but this choice is rather arbitrary. 

The first result to understand is a typical stress profile, without any friction fiuctuation, 
and the same parameters but with fiuctuation. This is shown in Figure 2, which displays the 
l/l/-component of the stress at the centerline of the bin and at the bin wall. For comparison, 
the Janssen stress is also shown. We have imposed the "initial condition" a = on y = 0. 
However, the condition for a surface y = h{x) to be stress free is (in general) inconsistent with 
y = const] imposing a = leads to a free boundary problem for the upper surface, a problem 
we do not wish to address here. The regular oscillations in Figure 2(a) are due to mismatch 
in the imposed stress at the intersection of the y = surface and the bin wall, and are well 
documented (see e.g. references in |]^); the period of these oscillations is related to the speed of 
the characteristics of the hyperbolic system. In Figure 2(b), fiuctuations at the walls are larger 
than at the centerline, and the wall stress is some 15% larger than centerline. Notice that the 
regular oscillations in Fig. 2(a) are dissipated by the randomness. 

In Figure 3, the yy-stiess is shown as a function of position across the bin, at the depth 
y = 10, the terminus of the computations in Figure 2. The variation across the bin illustrates 
the limitations of the Janssen assumptions. Nonetheless, Figures 2 and 3 show that the Janssen 
analysis provides a good estimate of the centerline stress (and NOT of the average stress!). This 
partially explains why the measurements of [^ are larger than the Janssen predictions. The 
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Figure 2: (a) The |/|/-component of the stress at the centerhne and the wall, with no random 
component of friction. For comparison the Janssen solution is also plotted. Here the nominal 
internal friction angle = 30° and the nominal wall friction angle 5 — 15°. (b) Similar plot, but 
with a random component added to the friction angles. 
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Figure 3: Variation in the yy-stiess across the width of the bin, at y = 10. Shown are results 
both with and without a random component of the friction angles. In both cases, the nominal 
friction angles are (f) — 30°, 6 — 15°. 
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Figure 4: (a) The yy-component of the stress at the centerhne for different nominal internal 
friction angles. Nominal wall friction is held fixed 6 — 15°. (b) The yy-component of the stress 
at the centerline for different nominal wall friction angles, with fixed nominal internal friction 
angle = 30°. 

centerline stress is typically 15 — 20% smaller than the largest stresses, found at the wall. 

Figure 4 illustrates the sensitivity of computations to changes in the nominal friction angles. 
In Fig. 4 (a), wall friction is held fixed while the nominal internal friction angle is varied 
from 15° to 30° (N.B. Recall that the random fiuctuation is 5% of the nominal angle). With 
lower internal friction, fluctuations become more pronounced. We conjecture that this is due to 
a lower friction angle transmitting a smaller fraction of stress (and of stress fluctuations) to the 
walls, leaving a larger fraction of stress (and of stress fluctuations) to be transmitted vertically. 
Notice too that, at the smallest friction angle, the regular oscillations of the stress reappear. 
When internal friction is held constant but wall friction is varied. Fig. 4 (b), the stress saturates 
deeper in the bin, and fluctuations are not apparent until after this saturation. We note that, 
with no wall friction, no weight is transferred to the bin walls and a hydrostatic stress results. 
Similarly, when periodic boundary conditions are imposed, a hydrostatic stress results. 

In Figure 5, fluctuations for two sets of friction angles are plotted. In each case, the equations 
were solved through y — 50. The centerline stress for 20 < y < 50 was extracted, and the average 
computed; this average should be the asymptotic value of the stress. A normalized deviation 
from the mean was found by subtracting the mean from the sample value, and dividing by the 
mean. For viewing, one signal is offset by 0.05. For the baseline case (p = 30°, S = 15°, the 
internal friction angle varies by about ±1.5° and the wall friction angle by about ±0.75°; the 
stress exhibits fluctuations of about ±4%. For the case = 30°, S — 5°, the internal friction 
angle again varies by about ±1.5° but the wall friction varies by only ±0.25°; the stress fluctuates 
about ±2.5%. 

Figure 6 provides a plot of spectral power for the baseline case (p = 30°, 6 = 15°. The stress 
was computed to a depth of y = 50; recall from Figure 2 that, for the given friction angles, the 
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Figure 5: Normalized fluctuations in the centerline yy-stiess for two pair of friction angles. Both 
signals are demeaned; for viewing, the top signal is vertically offset by 0.05. 

stress saturates well before y = 10. The centerline stress is sampled at every second timestep, 
from about y = 20 to y = 50. The power is computed using a Welch window with overlap, on 
the last 2560 sampled values. Shown is the (base 10 log of the) power for four variations: (i) 
no random component of friction; (ii) a random component of both internal and wall friction; 
[iii) a random component added to internal friction only; (iv) a random component added to 
wall friction only. The power for wall friction only lies atop the spectrum for wall and internal 
friction. The power for internal friction only deviates from these at lower wavenumber. Thus 
fluctuations in the stress are essentially due to a random component in the wall friction angle. 
From Figure 5, these fluctuations range up to about ±4% of the mean. Recall, this variation is 
based on about a 5% variation in the friction coefficient. The fluctuations reported in 0] are as 
large as 20%. This comparison suggests that a 15 — 20% variation in the friction coefficient is 
not an unreasonable parameter in stochastic models like the present. 

Analysis of the q-model shows that the number of occurrences of a fluctuation of a given 
size decays exponentially with size. Recent experiments 0| on short bins verify this flnding, for 
stresses larger than the mean; stresses smaller than the mean decay like a power law. Figure 7 
presents the distribution for the random friction model. The equilibrium equations were solved 
to a depth y = 50, and the centerline stress was recorded. Ten thousand realizations were made. 
The average over all realizations was calculated, subtracted from the sample value, and this 
difference was normalized by the average. Figure 7 is a histogram of these relative deviations. 
The distribution of fluctuations appears Gaussian, not exponential. 

4. Summary We have reexamined the Janssen analysis incorporating a random component 
of friction, solving for the mean and the second moment of the stress. For comparison, the 
nonlinear equilibrium equations for a Mohr-Coulomb material with random friction are solved 
numerically. The analysis suggests that fluctuations are signiflcant only after the stress begins 
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Figure 6: Log (base 10) of the spectral power for four variations of the base case. The variations 
are no random component of friction, a random component of both internal and wall friction, 
a random component added to internal friction only, and a random component added to wall 
friction only. The nominal = 30°, 5 = 15°. 
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Figure 7: A histogram of the normahzed deviation of the centerhne yy-stress from the mean, 
evaluated at y = 50. Friction parameters were = 30°, 6 = 15°. Computed for 10,000 
realizations, the distribution appears Gaussian, not exponential as predicted by the q-model. 
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to saturate, a finding consistent with the experimental work of 0. The primary contribution 
to stress fluctuations is randomness in the wall friction, a boundary condition. The fluctuations 
found in this model are set by a free parameter defining the magnitude of random friction. Our 
choice of this parameter results in fluctuations of about 5% of the mean stress, much less than 
the 15 — 20% found in experiments. 

The results presented here are in qualitative agreement with those of Socolar's a-model. His 
calculations incorporate stress balance in both horizontal and vertical directions, and a balance 
of angular momentum. The essential feature of the a-model is that particle friction transmits 
stress from particle to particle and, ultimately, to the walls of a bin. These stresses, and any 
stress fluctuation, are partially absorbed by the wall. In contrast, the q- model only considers 
vertical forces; stresses predicted by the q-model are more like hydrostatic forces, and there is 
no mechanism for dissipating fluctuations. 

A difficulty faced by all of these models is correlations. Experiments [Q] show chains of 
particles experiencing high stress (the frequency of which falls off exponentially with size). These 
pictures, and many other experiments, suggest that grain forces are correlated. However, we lack 
adequate information to introduce correlations into models in a meaningful way. Experimental 
results reported by measure static forces on short bins, and show no evidence for correlations. 
The question of whether there are correlations, and over what lengthscales are they important, is 
central to the entire formulation of a continuum framework for granular materials. Experimental 
and theoretical work is necessary to understand the nature of correlations. 

Mueth et al. also study the frequency of fluctuations in three dimensional systems. They 
find that, for fiuctuations larger than the mean, the frequency of fiuctuation of a given size 
decays exponentially with size of fiuctuations. For fiuctuations smaller than the mean, the 
decay follows a power law. Furthermore, their findings are largely unaffected by changes in the 
boundary friction. For purposes of comparison with this work, several factors are important. 
The experimental set-up has a depth-width aspect ratio of about 1 - 1.5. The glass beads and 
acrylic used in the experiment are very low friction materials, with both internal and wall friction 
angles about 10 — 15°. From the continuum perspective, stresses measured in this arrangement 
are hydrostatic-like. Walls do not support the bead pack, and even moderate changes in the wall 
friction would have only minor effects on stress measurements. We do not view these findings as 
invalidating the random friction model proposed here, at least not for engineering applications. 
The random packing model offers one possible explanation of these experimental findings. The 
mean stress for this model is given by Eqn. 3 in the limit a — 0, and equals pgy\ the standard 
deviation, Eqn. 6 in the a — limit, is epg^^y. For a short bin, ?/ ?a 1, and fiuctuations are 
on the order of e times the mean stress. Packing variations, interpreted as voids fraction, can 
range up to 20 — 30%. However even this model does not explain all the physics of small aspect 
ratio bins. 
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